Exact exchange optimized effective potential and self-compression of 

stabilized jellium clusters 
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In this work, we have used the exchange-only optimized effective potential in the self-consistent 
calculations of the density functional Kohn-Sham equations for simple metal clusters in stabilized 
jellium model with self-compression. The results for the closed-shell clusters of Al, Li, Na, K, and 
Cs with N =2, 8, 18, 20, 34, and 40 show that the clusters are 3% more compressed here than in 
the local spin density approximation. On the other hand, in the LSDA, neglecting the correlation 
results in a contraction by 1.4%. 
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I. INTRODUCTION 



The Kohn-Sham (KS)0 density functional theory 
(DFT)|2( is one of the most powerful techniques in elec- 
tronic structure calculations. However, the exact form 
of the exchange-correlation functional is still unknown 
and, in practice, one must use approximations. The ac- 
curacy of the predictions of the properties depends on 
how one approximates this functional. The simplest one 
is the local spin density approximation (LSDA) in which 
one uses the properties of the uniform electron gas lo- 
cally 0. This approximation is in principle appropriate 
for systems in which the variations of the spin densi- 
ties n a are sufficiently slow. For finite systems and sur- 
faces which are highly inhomogeneous, the generalized 
gradient approximation (GGA)[3j is more appropriate. 
In spite of the success of the LSDA and GGA, it is ob- 
served that in some cases these approximations fail to 
predict even qualitatively correct behaviors^ IE IE 0- 
On the other hand, appropriate self-interaction corrected 
versions of these approximations are observed to lead to 
correct behaviors pl|8|. These observations motivates one 
to use functionals in which the self-interaction contribu- 
tion is removed exactly. One of the functionals, which 
satisfies this constraint, is the exact exchange (EEX) or- 
bital dependent functional. Using the EEX functional 
leads to the correct asymptotic behavior of the KS po- 
tential as well as to correct results for the high density 
limit in which the exchange energy is dominated @ . Al- 
though neglecting the correlation effects in orbital depen- 
dent functionals fails to reproduce the dispersion forces 
such as the van der Waals forces ^E EH- the EEX in 
some respects is advanta geo us over the local and semi- 
local approximations [ill Il2fl. To obtain the local ex- 
change potential from the orbital dependent functional, 
one should solve the optimized effective potential (OEP) 
integral equation. Recently, Kummel and Perdew |Lflll4l ] 
have invented an iterative method which allows one to 
solve the OEP integral equation accurately and efficiently 
even for three dimensional systems. This method is used 
in this work. 



To simplify the cluster problem, one notes that the 
properties of alkali metals are dominantly determined 
by the delocalized valence electrons. In these met- 
als, the Fermi wavelengths of the valence electrons are 
much larger than the metal lattice constants and the 
pseudo-potentials of the ions do not significantly affect 
the electronic structure. This fact allows one to replace 
the discrete ionic structure by a homogeneous positive 
charge background which is called jellium model (JM). 
In its simplest form, one applies the JM to metal clus- 
ters by replacing the ions of an iV-atom cluster with 
a sphere of uniform positive charge density and radius 
R = (zA r ) 1 / 3 r s , where z is the valence of the atom and 
r s is the bulk value of the Wigner-Seitz (WS) radius for 
valence electrons ^E EE El- Assuming the spherical ge- 
ometry is justified only for closed-shell clusters which is 
the subject in this work. However . it is a known fact 
that the JM has some drawbacks (lit Il9| . The stabilized 
jellium model (SJM) in its original form|20| was the first 
attempt to overcome the deficiencies of the JM and still 
keeping the simplicity of the JM. Application of the SJM 
to simple metals and metal clusters has shown signif- 
icant improvements over the JM results 20]. However, 
for small metal clusters the surface effects are important 
and the cluster is self-compressed due to its surface ten- 
sion. This effect has been successfully taken into account 
by the SJM which is called SJM with self-compression 
(SJM-SC)[2l],|22. Application of the LSDA-SJM-SC to 
neutral metal clusters has shown that the equilibrium 
r s values of small clusters are smaller than their bulk 
counterparts and approaches to it for very large clusters. 
This trend is consistent with the results of ab. initio. 
calculations^!!^. 

In this work we have used the EEX-SJM-SC to ob- 
tain the equilibrium sizes and energies of closed-shell 
neutral A^-electron clusters of Al, Li, Na, K, and Cs 
for N=2, 8, 18, 20, 34, and 40 (for Al, N = 18 cor- 
responds to Alg cluster and other values do not corre- 
spond to a real Al n ). In order to have an estimate for 
the self-interaction effects, we have repeated the calcu- 
lations for exchange-only local spin density approxima- 
tion (x-LSDA) in which the spin-polarized version of the 



2 



Dirac form, E x = c x J dr n 4 / 3 , is used. Comparison of 
the results shows that (except for N = 40 in Al case) 
the relation f EEX < f x -LSDA < r hSDA - The organiza- 
tion of this paper is as follows. In section [H] we explain 
the calculational schemes. Section Hill is devoted to the 
results of our calculations and finally, we conclude this 
work in section llVl 



II. CALCULATIONAL SCHEMES 

In this section we first explain how to implement the 
exact exchange in the SJM, and then will explain the 
procedure for the OEP calculations. 

A. Exact exchange stabilized jellium model 
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As in the original SJM 20] , here the Ashcroft empty 
core pseudo-potential[25| is used for the interaction of an 
ion of charge z with an electron at a relative distance r: 



, , ( -2z/r , (r > r c ) 
W{r) = { , \r<r') 



(1) 



The core radius, r c , will be fixed by setting the pres- 
sure of the bulk system equal to zero. In the EEX-SJM, 
the average energy per valence electron in the bulk with 
density n is given by 



e(n) = t s (n) + e x (n) +w R (n,r c ) +e M (n), (2) 



with 



t s (n) = c k n 



£ x (n) 



c x n 



1/3 



Ck 



:(3tt 



2\2/3 



:(3A 



(3) 



(4) 



(5) 



All equations throughout this paper are expressed in 
Rydberg atomic units. Here t a and e x are the kinetic 
and exchange energy per particle, respectively, wr is 
the average value of the repulsive part of the pseudo- 
potential (wr = Aimr^), and £m is the average Madelung 
energy. Demanding zero pressure for the bulk system at 
equilibrium yields: 



0. (6) 



{2t s (n) + e x {n) + Yl-nnr\ + e M {n)} n=n 
Solution of this equation for r c gives 

r c (rf ) = {-2t s (r s ) - e x {r s ) ~ e M (r s )}^ r 



(7) 



FIG. 1: Pseudo-potential core radii in atomic units for dif- 
ferent rf values. 



In Fig. [I]we have plotted the core radii for different 
values of rf which assume 2.07, 3.28, 3.99, 4.96, and 
5.63 for Al, Li, Na, K, and Cs, respectively. The result is 
compared with the case in which the correlation energy 
is also incorporated (see Eq.(26) of Ref.[20|). As is seen, 
to stabilize the bulk system in the EEX case, the core 
radii assume smaller values. 

As in the original SJM 20] (but in the absence of the 
correlation energy component), at equilibrium density we 
have 



(Sv) 



1 



ws 



[2t s {n B ) + e x (n B )}. 



(8) 



Here, {8v)ws is the average of the difference potential 
over the WS cell and the difference potential, Sv, is de- 
fined as the difference between the pseudo-potential of a 
lattice of ions and the electrostatic potential of the jel- 
lium positive background. Once the values of (Sv)ws 
and r c as functions of rf are found, the EEX-SJM total 
energy of a cluster becomes 

EEEXSJM[n-[,ni,r s ,rf] = E E EX-JM[n-\,ni,r s ] 



+(e M + w R ) J dr n+ (r) 
+ (fo)ws I <fre(r)[n(r)-n+(r)]. (9) 



Here. 



EEEX-jM[n h n u r s ] = T s [n T ,?iJ + E x [nf,n{\ 

^ dr 0([n,n+];r) [n(r) -n+(r)], (10) 



dr dr 



(11) 
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^n,n + ];r) = 2l d r ^ r '\- n :P , (12) 



r r 



i W= E EiM') 

<r=T4 i=l 



?i + (r) = n6(R-r); 
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(13) 



(14) 



To obtain the equilibrium size and energy of an N- 
atom cluster in EEX-SJM-SC, we solve the equation 



^-E(N,r Sl r c ) 



= 0, 



(15) 



r s =f s (N) 



where N and r c are kept constant and E is given by Eq. 
(jSJ . The preocedure for the x-LSDA is the same as above 
except for that the Dirac exchange energy must be used. 



B. The OEP equations 

Kiimmel and Perdew[l4) have proved, in a simple way, 
that the OEP integral equation is equivalent to 



^V4(r)<Mr)+c.c. = 0. 



(16) 



i=l 



ipiu are the self-consistent KS orbitals and ipi a are orbital 
shifts. The self-consistent orbital shifts and the local ex- 
change potentials are obtained from the iterative solu- 
tions of inhomogeneous KS equations. Taking spherical 
geometry for the jellium background and inserting 



and 



Mr) = *^Y htmt ((l), 



Mr) = ^2lWfi), 



(17) 



(18) 



in to the inhomogeneous KS equation (Eq.(21) of 
Ref.Q) one obtains 26] 



-p +e ia -v effa (r) 



k{li + 1) 



Here, e% a are the KS eigenvalues and 



v e ff a (r) = v(r) + v H (r) + v xa (r), 



vh (r) = 2 J dr 



n(r') 



(19) 



(20) 



(21) 



The right hand side of Eq. l|19f) can be written as 



qiv{r)=qt\r)+q\;\r), 



(22) 



with 



1il\ r ) = [Vxca(r) ~ Vxda + U xcia ] Xia(r), (23) 



and 



j=l l=\U~li\ 

x[I(ljrrij,limi,lmj — mi)] 2 . (24) 
The quantities B and / in Eq. I|24(l are defined as 
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r 

r T+i 



+ / dr'xiair')^')^ (25) 



imm^kmulm) = / dtt Y t * (n)Y hmi (Q.)Y lm (Q,), 



(26) 

and the bar over I 2 implies average over m, and rrij. 
Also, the expression for u X i a reduces to 

4tt 



n„ U+lj 

z «- = - 2 E E 2/TT 

j=l l=\U-l 3 \ 



[I(ljmj,kmi, Irrij - mi)]" 



x / dr Xi<T{r)Xi<T{r)B„{i,j,l,r). (27) 
Jo 

The procedure for the self-consistent iterative solutions 
of the OEP equations is explained in Refs.0,|2(j. 

In Fig. |21 the self-consistent source terms qi a (r) of 
Eq. (|19fl are plotted for the equilibrium size of Nais 
cluster. The corresponding orbital shifts (r) are shown 
in Fig|3 



III. RESULTS AND DISCUSSION 

We have used the EEX-SJM-SC to obtain the equilib- 
rium sizes and energies of closed-shell 2, 8, 18, 20, 34, 
and 40-electron neutral clusters of Al, Li, Na, K, and Cs. 

In Tabled we have listed the equilibrium r s values, to- 
tal energies and exchange energies. As is seen, the equi- 
librium r s values of the clusters are almost the same up to 
3 decimals for the KLI and OEP schemes whereas, there 
are significant differences between the OEP, x-LSDA, and 
LSDA values. As an example, we have plotted the equi- 
librium r s values of the closed-shell Kjv clusters in Fig. 
0] It shows that the LSDA predicts larger cluster sizes 
than the x-LSDA and OEP. 
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FIG. 2: Right hand side of Eq. I|19p for the self-consistent 
equilibrium size of Naig. 



FIG. 4: Equilibrium r a values of Kjv clusters for different 
sizes. The dotted line is rf. KLI and OEP predict smaller 
sizes for the clusters than the x-LSDA and LSDA. 
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FIG. 3: Orbital shifts in atomic units for the self-consistent 
equilibrium size of Naig. 



To illustrate the trend in the r s values, we plot the 



difference ( 



tLSD 



xKLI 



) for all species in Fig. |SJ One 



notes that for a given element, the difference is larger for 
smaller clusters. On the other hand, the difference for the 
lower-density element is higher. However, the difference 
is about 3% on average. We therefore conclude that the 
EEX-SJM-SC predicts smaller bond lengths compared to 
the LSDA-SJM-SC. Comparison of the f s values for the 
LSDA and x-LSDA shows that bond lengths in the LSDA 
is about 1.4% larger on average. This difference should be 
attributed to the correlation effects. On the other hand, 
the same comparison between x-LSDA and KLI shows 
that, except for N = 40 in Al, f x ^LSDA > ^kli by 1.5% 
on average. This difference is due to the self-interaction 
effects in the Dirac form for the exchange functional. 




Q 
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N 

FIG. 5: Difference in the equilibrium r s values of clusters 
with different sizes. The difference is larger for lower-density 
elements. 



Comparison of the equilibrium total energies of the 
OEP and KLI shows that OEP energies are on average 
0.02% more negative. This result should be compared to 
the simple JM results pf3| which is 1.2%. On the other 
hand, comparison of the exchange energies shows that on 
the average, the exchange energies in OEP is 0.34% more 
negative than those in the KLI. 

In Table [HJ we have listed the lowest and highest oc- 
cupied KS eigenvalues for different schemes. As in the 
simple JM [26|. the OEP KS eigenvalue bands are con- 
tracted relative to those of the KLI. That is, for all N, the 
relation A OEP < A KLI holds. Here, A = e H - £l is the 
difference between the maximum occupied and minimum 
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occupied KS eigenvalues. For the same external poten- 
tial, the OEP and KLI results coincide for two-electron 
systems and A = 0. The results in TablJTflshow that the 
maximum relative contraction, |A OEP — A KLI \/ A KLI , is 
2.7% which corresponds to Csig- 

The same comparisons between OEP and x-LSDA 
shows that E OEP < E X ~ LSDA by 5.2% on average, and 
E° EP < E%- LSDA by 11% on average. The band widths 
do not show any regular pattern, however, in the OEP 
the bands mostly contract relative to the x-LSDA. 

Finally, we compare the results of LSDA and x-LSDA, 
which will show the correlation effects. As is seen in 
Table U the total energies are close to each other for 
the high-density cases. That is, in the high density limit 
the exchange dominates the correlation. However, the 
total energies in the LSDA are more negative by 10% 
on average which is due to the correlation effects. On 
the other hand, the difference in the exchange energies is 
about 0.96% on average which is quite a small fraction. In 
the high density limit, the inequality E*~ LSD < E LSDA 
holds whereas, in the low density limit the inequality 
changes sign. 

IV. SUMMARY AND CONCLUSION 

In this work, we have considered the exact-exchange 
stabilized jellium model with self-compression in which 
we have used the exact orbital-dependent exchange func- 



tional. This model is applied for the simple metal clusters 
of Al, Li, Na, K, and Cs. For the local exchange poten- 
tial in the KS equation, we have solved the OEP integral 
equation by the iterative method. By finding the min- 
imum energy of an A-atom cluster as a function of r s , 
we have obtained the equilibrium sizes and energies of 
the closed-shell clusters (N = 2,8,18,20,34,40) for the 
four schemes of LSDA, KLI, OEP, and x-LSDA. The re- 
sults show that in the EEX-SJM, the clusters are more 
contracted relative to the x-LSDA-SJM, i.e., 1.5% more 
contraction on average. The KLI and OEP results show 
equal values (up to three decimals) for the equilibrium 
r s values. The equiliblium sizes in LSDA and x-LSDA 
differ by 1.4% on average. In the LSDA and KLI the dif- 
ference in 3% on average. The total energies in the OEP 
are more negative than the KLI by 0.02% on the average. 
It should be mentioned that in the simple JM the KLI 
and OEP total energies for Al were positive (except for 
N = 2). On the other hand, the exchange energies in the 
OEP is about 0.34% more negative than that in the KLI. 
Comparison of the OEP and x-LSDA shows a duifference 
of 5.2% in the total energies and 11% in the exchange. 
The difference in the exchange energies of LSDA and x- 
LSDA is small (about 0.96%) whereas the total energy 
in the LSDA is about 10% more negative which is due to 
the correlation effects. The widths of the occupied bands, 
Eh — £l in the OEP are contracted relative to those in 
the KLI by at most 2.7%. 
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TABLE I: Equilibrium sizes, r s , in bohrs, the absolute values of total and exchange energies in rydbergs are compared for KLI, 
OEP,x-LSDA, and LSDA schemes. In the LSDA, the total energies include the correlation energies as well. 
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Al" 2.07 
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1.430 


1.5700 


0.9253 
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1.4364 


0.7574 
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7 


.4223 


40 


3 


.813 
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.0432 


3.813 


16.1431 


9 


.0632 


3.843 
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3. 


.855 
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.3920 


4. 


354 
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3 


.3860 


4.781 


5.8178 


3. 


1269 


4. 


805 


6.6130 


3 


1427 


20 


4. 


.710 


6.6722 


3 


.6887 


4. 


710 


6.6733 


3 


.7002 


4.760 


6.4820 


3. 


.4670 


4. 


784 


7.3628 


3 


.4795 


34 


4. 


787 


11.3534 


6 


.3392 


4. 


787 


11.3579 


6 


.3782 


4.824 


11.0656 


5 


.9836 


4. 


.840 


12.5823 


6 


.0079 


40 


4. 


.768 


13.2650 


7 


.2975 


4. 


768 


13.2671 


7. 


.3162 


4.808 


13.0090 


7 


.0029 


4. 


.824 


14.7863 


7 


.0226 



Cs 5.63 



2 


5.006 


0.6123 





.3494 


5.006 


0.6123 


0.3494 


5.109 


0.5624 





.2856 


5.213 


0.6395 





.2910 


8 


5.261 


2.3990 


1 


.3322 


5.261 


2.3994 


1.3363 


5.342 


2.2918 


1. 


2039 


5.383 


2.6135 


1 


.2133 


18 


5.390 


5.3547 


2 


.9775 


5.389 


5.3564 


2.9963 


5.443 


5.1842 


2. 


.7658 


5.472 


5.9215 


2 


.7821 


20 


5.366 


5.9403 


3 


.2589 


5.366 


5.9414 


3.2701 


5.425 


5.7729 


3. 


0640 


5.451 


6.5894 


3 


.0784 


34 


5.445 


10.1156 


5 


.6044 


5.445 


10.1200 


5.6423 


5.488 


9.8599 


5. 


.2875 


5.508 


11.2652 


5 


.3122 


40 


5.428 


11.8123 


6 


.4416 


5.428 


11.8144 


6.4598 


5.472 


11.5881 


6 


.1873 


5.493 


13.2347 


6 


.2061 



"Here, N=18 corresponds to Al6 cluster and other Af's do not 
correspond to a real Al clusters. 



7 



TABLE II: The absolute values at equilibrium state of the highest occupied and lowest occupied Kohn-Sham eigenvalues in 
rydbergs are compared for KLI, OEP, x-LSDA, and LSDA schemes. 









KLI 




OEP 




-LSDA 




LSDA 


Atom 


N 


£ z, 


£ H 


el 


EH 


£ L 




£l 


£i? 


Al 


2 


8152 


8152 


8152 


8152 


n a ^ f\7 

U.^iOO / 


4367 


5012 


5012 




g 


1 1142 


6714 


1 1088 


6713 


8201 


3919 


8821 


4605 




18 


1.1727 


0.5507 


1.1619 


0.5492 


0.9497 


0.3310 


1.0129 


0.4009 




20 


1.1856 


0.5000 


1.1804 


0.4993 


0.9665 


0.2964 


1.0282 


0.3622 




34 


1.2055 


0.4826 


1. 1998 


0.4789 


1.0192 


0.2939 


1.0853 


0.3649 




40 


1 2202 


D 44Q0 

U . tit: f \-J 


1 2136 


4450 


1 0541 


2761 


1 0965 


3401 


Li 


2 


4777 


4777 


4777 


4777 


2445 


2445 


2983 


2983 




g 


5760 


4157 


5735 


4158 


n ^q^^ 
u . oyoo 


2383 


4476 


2937 




18 


0.5935 


0.3601 


0.5879 


0.3591 


0.4523 


0.2179 


0.5062 


0.2738 




20 


0.5889 


0.3221 


0.5865 


0.3228 


0.4522 


0.1893 


0.5061 


0.2423 




34 


6029 




5991 


n 39^ a 
u . ozoo 


4832 


n 90^ 


5374 


O 9fi9D 




40 


5979 


2971 


u . oyou 


2960 


4816 


1832 


u.oooo 


U.ZOOD 


In a 


2 


yj.oooo 




y.).oooo 


u.oooo 




u. iyoi 


n 9/1 T7 


fl 9/1 




g 


4467 


3406 


4451 


3408 


u . zy do 


1936 


3453 


2434 




18 


0.4544 


0.2989 


0.4502 


0.2981 


0.3373 


0.1805 


0.3859 


0.2308 




20 


0.4485 


0.2672 


0.4470 


0.2682 


0.3361 


0.1565 


0.3851 


0.2042 




34 


4583 


2750 


4551 


2730 


U.OOOO 


1728 


4078 






40 


4520 


2480 


4502 


2477 


3564 


1534 


4055 


2017 


K 


2 


3100 


3100 


3100 


3100 


1526 


1526 


1957 


1957 




g 


3408 


2733 


u . ooyo 


2735 


2188 


1536 


n 9R99 


1977 




18 


0.3416 


0.2422 


0.3385 


0.2415 


0.2463 


0.1458 


0.2894 


0.1902 




20 


0.3356 


0.2169 


0.3349 


0.2182 


0.2450 


0.1266 


0.2885 


0.1687 




34 


0.3419 


0.2246 


0.3393 


0.2230 


0.2607 


0.1413 


0.3042 


0.1861 




40 


0.3355 


0.2025 


0.3346 


0.2027 


0.2585 


0.1255 


0.3023 


0.1683 


Cs 


2 


0.2723 


0.2723 


0.2723 


0.2723 


0.1324 


0.1324 


0.1724 


0.1724 




8 


0.2923 


0.2405 


0.2913 


0.2406 


0.1843 


0.1343 


0.2247 


0.1752 




18 


0.2907 


0.2141 


0.2880 


0.2134 


0.2061 


0.1285 


0.2461 


0.1696 




20 


0.2850 


0.1921 


0.2847 


0.1935 


0.2048 


0.1118 


0.2454 


0.1509 




34 


0.2897 


0.1992 


0.2873 


0.1977 


0.2176 


0.1252 


0.2580 


0.1668 




40 


0.2835 


0.1797 


0.2831 


0.1801 


0.2157 


0.1115 


0.2564 


0.1512 



